######
# This code reproduces the analysis in
# Ahlquist/Clayton/Levi
# International Organization 2014
# LAST UPDATED: 19 Aug 2013
# UPDATED BY: ABC
# NOTE: this uses Amelia 1.6.3, ZeligChoice 0.7-0 and Zelig 4.1-2
######

rm(list=ls())
##load libraries 
library(MatchIt)
library(car)
library(vcd)
library(VGAM)
library(MASS)
library(Amelia)
library(foreign)
library(modeest)
library(ineq)


setwd("/Users/amandaclayton/Desktop/DataVerse")
modeldata<-read.csv("ACLTradeAttitudes.csv")
modeldata<-modeldata[,2:21]
names(modeldata)

##Use Amelia to impute missing data
bds1 <- matrix(c(4, 18, 96), nrow = 1, ncol = 3)  #set bounds for age (actual range in data)
bds2 <- matrix(c(2, 0, 52), nrow = 1, ncol = 3)  #set bounds for ilwu tenure (actual range in data)
bds3<-rbind(bds1, bds2)
 
m<-20
a.out <- amelia(modeldata, m = m,  ord=c("income", "edu", "round","nafta", "ilwu_fam"), nom=c("mode", "local", "imports_cat", "pol_party","white", "female","A", "B", "ilwu",  "hispanic", "union_fam", "union"), bounds=bds3)

##average across all 20 datasets
## no missingness in ILWU and Local, date
## mode for white, female, union family member and hispanic
## mean for age, education, income, nafta, ilwu tenure  

#create age average

temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$age}
comb.age<-temp/m

#create ilwu tenure average
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$ilwu_years}
comb.ilwu_yrs<-round(temp/m,1)

#create education average
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$edu}
comb.edu<-round(temp/m,0)

#create income average 
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$income}
comb.income<-round(temp/m,0)

#create support for NAFTA average 
temp<-0
for(i in 1:m){
  temp<-temp+a.out$imputations[[i]]$nafta}
comb.nafta<-round(temp/m,0)

#pol party mode (used for specification checks)
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$pol_party)}
party.mode<-apply(temp, 1, mfv)

#randomly filling in value where there are multiple modes
for(i in 1:length(party.mode)){
  if(sum(party.mode[[i]]==c(1,2))==2) party.mode[[i]]<-sample(c(1,2),1)
  if(sum(party.mode[[i]]==c(2,1))==2) party.mode[[i]]<-sample(c(1,2),1)
  if(sum(party.mode[[i]]==c(2,3))==2) party.mode[[i]]<-sample(c(3,2),1)
  if(sum(party.mode[[i]]==c(3,2))==2) party.mode[[i]]<-sample(c(3,2),1)
  if(sum(party.mode[[i]]==c(3,1))==2) party.mode[[i]]<-sample(c(3,1),1)
  if(sum(party.mode[[i]]==c(1,3))==2) party.mode[[i]]<-sample(c(3,1),1)
}
party.mode<-as.numeric(party.mode)

##union family member
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$union_fam)}
union_fam.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(union_fam.mode)){
  if(sum(union_fam.mode[[i]]==c(0,1))==2) union_fam.mode[[i]]<-sample(c(0,1),1)
  if(sum(union_fam.mode[[i]]==c(1,0))==2) union_fam.mode[[i]]<-sample(c(0,1),1)
}
union_fam.mode<-as.numeric(union_fam.mode)

##white
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$white)}
white.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(white.mode)){
  if(sum(white.mode[[i]]==c(0,1))==2) white.mode[[i]]<-sample(c(0,1),1)
  if(sum(white.mode[[i]]==c(1,0))==2) white.mode[[i]]<-sample(c(0,1),1)
}
white.mode<-as.numeric(white.mode)

###female
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$female)}
female.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(female.mode)){
  if(sum(female.mode[[i]]==c(0,1))==2) female.mode[[i]]<-sample(c(0,1),1)
  if(sum(female.mode[[i]]==c(1,0))==2) female.mode[[i]]<-sample(c(0,1),1)
}
female.mode<-as.numeric(female.mode)

##union member
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$union)}
union.mode<-apply(temp, 1, mfv)

#randomly filling in value where there are multiple modes
for(i in 1:length(union.mode)){
  if(sum(union.mode[[i]]==c(0,1))==2) union.mode[[i]]<-sample(c(0,1),1)
  if(sum(union.mode[[i]]==c(1,0))==2) union.mode[[i]]<-sample(c(0,1),1)
}
union.mode<-as.numeric(union.mode)



##DV imports cat
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$imports_cat)}
imports.mode<-apply(temp, 1, mfv)

for(i in 1:length(imports.mode)){
  if(sum(imports.mode[[i]]==c(1,2))==2) imports.mode[[i]]<-sample(c(1,2),1)
  if(sum(imports.mode[[i]]==c(2,1))==2) imports.mode[[i]]<-sample(c(1,2),1)
  if(sum(imports.mode[[i]]==c(2,3))==2) imports.mode[[i]]<-sample(c(3,2),1)
  if(sum(imports.mode[[i]]==c(3,2))==2) imports.mode[[i]]<-sample(c(3,2),1)
  if(sum(imports.mode[[i]]==c(3,1))==2) imports.mode[[i]]<-sample(c(3,1),1)
  if(sum(imports.mode[[i]]==c(1,3))==2) imports.mode[[i]]<-sample(c(3,1),1)
  if(sum(imports.mode[[i]]==c(0,1))==2) imports.mode[[i]]<-sample(c(0,1),1)
  if(sum(imports.mode[[i]]==c(1,0))==2) imports.mode[[i]]<-sample(c(0,1),1)
  if(sum(imports.mode[[i]]==c(0,2))==2) imports.mode[[i]]<-sample(c(0,2),1)
  if(sum(imports.mode[[i]]==c(2,0))==2) imports.mode[[i]]<-sample(c(0,2),1)
  if(sum(imports.mode[[i]]==c(0,3))==2) imports.mode[[i]]<-sample(c(0,3),1)
  if(sum(imports.mode[[i]]==c(3,0))==2) imports.mode[[i]]<-sample(c(0,3),1)
}

imports.mode<-as.numeric(imports.mode)

###hispanic
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$hispanic)}
hispanic.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(hispanic.mode)){
  if(sum(hispanic.mode[[i]]==c(0,1))==2) hispanic.mode[[i]]<-sample(c(0,1),1)
  if(sum(hispanic.mode[[i]]==c(1,0))==2) hispanic.mode[[i]]<-sample(c(0,1),1)
}
hispanic.mode<-as.numeric(hispanic.mode)

###A-status
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$A)}
A.mode<-apply(temp, 1, mfv)

#ranomly filling in value where there are multiple modes
for(i in 1:length(A.mode)){
  if(sum(A.mode[[i]]==c(0,1))==2) A.mode[[i]]<-sample(c(0,1),1)
  if(sum(A.mode[[i]]==c(1,0))==2) A.mode[[i]]<-sample(c(0,1),1)
}
A.mode<-as.numeric(A.mode)

###B-status
temp<-NULL
for(i in 1:m){
  temp<-cbind(temp,a.out$imputations[[i]]$B)}
B.mode<-apply(temp, 1, mfv)

#randomly filling in value where there are multiple modes
for(i in 1:length(B.mode)){
  if(sum(B.mode[[i]]==c(0,1))==2) B.mode[[i]]<-sample(c(0,1),1)
  if(sum(B.mode[[i]]==c(1,0))==2) B.mode[[i]]<-sample(c(0,1),1)
}
B.mode<-as.numeric(B.mode)

comb.data<-as.data.frame(cbind(modeldata$ilwu, modeldata$local, round(a.out$imputations[[1]]$date, 0), round(comb.age, 0), white.mode, comb.edu, female.mode, hispanic.mode, A.mode, B.mode, party.mode, comb.income, a.out$imputations[[1]]$mode, union.mode, comb.ilwu_yrs, union_fam.mode, imports.mode , comb.nafta))

names(comb.data)<-c("ilwu", "local", "date", "age", "white", "edu", "female", "hispanic", "A", "B", "pol_party", "income", "mode", "union", "ilwu_yrs", "union_fam", "imports_cat", "nafta" )

#recode factor variables as factors
comb.data$imports_cat<-as.factor(comb.data$imports_cat)
comb.data$pol_party<-as.factor(comb.data$pol_party)
comb.data$local<-as.factor(comb.data$local)
comb.data$mode<-as.factor(comb.data$mode)

###Matching
m.out1<-matchit(ilwu~white+as.factor(local)+age+edu+female+ date+hispanic,
                method="genetic", discard="both",data = comb.data, ratio=2)
match.out<-match.data(m.out1)
##matching diagnostics
summary(m.out1)

#This code produces the plots in the matching appendix

#This code produces Figure A1
##100 x the mean diff/pooled SD 
summ.plot <- summary(m.out1, standardize=TRUE) 
plot(summ.plot, interactive=F, lines=NULL)
text(x=rep(0.9, 9), y=c(1.95, 1.55, 0.86, 0.74, 0.68, 0.35, 0.24, 0.158, 0.115, 0.065),
     labels=c("distance", "female", "age", "education", "Seattle", "Tacoma", "Long Beach", "Hispanic", "white", "date"), cex=0.7)


#This code produces Figure A2
plot(m.out1, type="hist") 

#This code produces Figure A3
plot(m.out1, type="jitter", interactive=F) 
summary(match.data(m.out1))


############# Trade opinion models for both IO paper & Ahlquist/Levi 2013 ###########

####This code produces Figure 2
##split data ILWU/RDD
ilwudata<-subset(modeldata4, ilwu==1) #638
rdddata<-subset(modeldata4, ilwu==0) #604
ilwu_unmatched<-table(ilwudata$imports_cat2)/length(na.omit(ilwudata$imports_cat2))
 #not imputed 
rdd_unmatched<-table(rdddata$imports_cat2)/length(na.omit(rdddata$imports_cat2))
 #not imputed 
##to compare w/ ANES data (from the 2004 ANES time series study)

nes<-read.dta("anes2004TS.dta")
imports<-nes$V045114
imports<-as.character(imports)
imports[imports=="1. Favor"]<-1 
imports[imports=="5. Oppose"]<-2 
imports[imports=="7. Haven't thought much about this"]<-3 
imports[imports=="8. Don't know"]<-NA 
imports[imports=="9. Refused"]<-NA
ANES<-table(imports)/(length(na.omit(imports)))

unmatched_table<-(rbind(ilwu_unmatched, rdd_unmatched, ANES))

#for matched data
##split matched data ILWU/RDD
MatchILWUdata<-subset(match.out, ilwu==1) #637
MatchRDDdata<-subset(match.out, ilwu==0) #256
ilwu_matched<-table(MatchILWUdata$imports_cat)/length(na.omit(MatchILWUdata$imports_cat))
rdd_matched<-table(MatchRDDdata$imports_cat)/length(na.omit(MatchRDDdata$imports_cat))

matched_table<-rbind(ilwu_matched, rdd_matched)

par(mfrow=c(1,2))
barplot(unmatched_table, beside=T, main="Support for Import Restrictions \n Unmatched Sample", ylim=c(0, 0.6),
        xlab=c(""), xaxt="n", col=c("black", "dark grey", "white"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1))
axis(side=1, at=c(0.35, 1.2, 2.2), labels=c("Favor", "Oppose", "Haven't Thought \n About It"),
     tick=FALSE, cex.axis=0.7)
legend("topleft", fill= c("black", "dark grey", "white"), border="black", col=c("black", "dark grey", "white"),
       legend=c("ILWU", "RDD", "ANES"), cex=0.6)

barplot(matched_table, beside=T, main="Support for Import Restrictions \n Matched Sample", ylim=c(0, 0.6), xlab=c(""),
        xaxt="n", col=c("black", "dark grey"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1))
axis(side=1, at=c(0.35, 1, 1.6), labels=c("Favor", "Oppose", "Haven't Thought \n About It"), tick=FALSE, cex.axis=0.7)
legend("topleft", fill=c("black", "dark grey"), border="black", col=c("black", "dark grey"),
       legend=c("ILWU", "RDD"), cex=0.6)

###This code produces Table 2 (IO)
##main imports model
##NOTE: income, PID, survey mode all excluded as post-treatment variables
#
impres.1<-vglm(as.factor(imports_cat)~ ilwu + A + B +
              white + female + as.factor(local) +
              age + I(age^2) + edu + date + hispanic,
      family = multinomial,data=match.out, weights = match.out$weights) 
summary(impres.1)
yhats<-predictvglm(impres.1, type="response")
#with alternate reference category:
alt<-as.data.frame(match.out)
##original coding: 1=favor, 2=oppose, 3=dk (ref cat)
alt$imports_cat<-recode(alt$imports_cat, "1=2; 2=3; 3=1") #1=don't know, 2=favor, 3=oppose (ref cat)
impres.alt<-vglm(as.factor(imports_cat)~ ilwu + A + B +
              white + female + as.factor(local) +
              age + I(age^2) + edu + date + hispanic,
      family = multinomial,data=alt, weights = alt$weights) 
summary(impres.alt)

##Simulations and ternary plots
#calculating risk ratio 
x.trt.A<- c(1,1,0,1,0,3,mean(match.out$age), mean(match.out$edu),
            round(mean(match.out$date),0),0)
x.cntrl<- c(0,0,0,1,0,3,mean(match.out$age), mean(match.out$edu),
            round(mean(match.out$date),0),0)
x.trt.A<- data.frame(t(x.trt.A))
x.cntrl<- data.frame(t(x.cntrl))

names(x.trt.A)<-names(x.cntrl)<-c("ilwu", "A", "B", "white", "female", "local",
              "age", "edu" , "date" , "hispanic")
risk.rat<-predictvglm(impres.1, type="response", newdata=x.trt.A)/predictvglm(impres.1, type="response", newdata=x.cntrl)


y.hat.rdd<-yhats[match.data(m.out1)$ilwu==0,]
y.hat.ilwu<-yhats[match.data(m.out1)$ilwu==1,]
y.hat.A<-yhats[match.data(m.out1)$A==1,]
y.hat.B<-yhats[match.data(m.out1)$B==1,]
y.hat.C<-yhats[match.data(m.out1)$A==0 & match.data(m.out1)$B==0 & match.data(m.out1)$ilwu==1,]
y.hats.ABC<-rbind(y.hat.A,y.hat.B, y.hat.C)
A.men<-c(rep(1,nrow(y.hat.A)), rep(0,nrow(y.hat.B)), rep(0,nrow(y.hat.C)))
B.men<-c(rep(0,nrow(y.hat.A)), rep(1,nrow(y.hat.B)), rep(0,nrow(y.hat.C)))
y.hats.ABC<-cbind(y.hats.ABC,A.men, B.men)

##This code produces Figure 3
ternaryplot(
  yhats,
  pch = ifelse(match.out$ilwu == 0, 16, 2),
  col = ifelse(match.out$ilwu == 0, "black" , grey(0.5)),
  main="Predicted Probabilities \n ILWU/RDD Attitudes on Placing Limits on Imports",
  dimnames=c("favor","oppose","haven't thought about it"),
  cex = .7
)
grid_legend(0.2,0.63, pch=c(2,16), col=c(grey(0.5), "black"), labels=c("ILWU observations", "RDD observations"),
            title="", vgap=1, frame=F)

##This code produces Figure 4
###A men vs casuals

###A men vs casuals
ternaryplot(
   y.hats.ABC[,1:3],
  pch = ifelse(A.men==1 & B.men!=1, 17,
  		ifelse(A.men!=1 & B.men==1, 3, 
  		1)),
  col = ifelse(A.men==1 & B.men!=1, grey(0.7),
  		ifelse(A.men!=1 & B.men==1, grey(0.5), 
  		"black")),
  main="Predicted Probabilities \n A/B/Casuals' Attitudes on Placing Limits on Imports",
  dimnames=c("favor","oppose","haven't thought about it"),
  cex = .6
)


grid_legend(0.2,0.63, pch=c(17,3,1), col=c(grey(0.7), grey(0.5), "black"), labels=c("Class-A", "Class-B", "Casuals"),
            title="", vgap=1, frame=F)


#### compare NAFTA responses to WSV data
#this code produces Figure 6
setwd("/Users/amandaclayton/Documents/labor/trade attitudes/other surveys")
wvs<-read.dta("wvs.dta") ##from 2006, most recent tabulated year of WVS
wvsUS<-as.data.frame(subset(wvs, wvs$s024a=="united states (5)"))

nafta<-wvsUS$e069_24
table(nafta)
nafta1<-recode(nafta, "'a great deal' = 1; 'quite a lot' = 2; 'not very much' = 3; 'none at all' = 4; NA=NA")
wvs_nafta<-table(nafta1)/length(na.omit(nafta1))  #n = 1170
ilwu_nafta<-table(ilwudata$nafta)/length(na.omit(ilwudata$nafta))
 #not imputed 
rdd_nafta<-table(rdddata$nafta)/length(na.omit(rdddata$nafta))
 #not imputed 

naftatable<-rbind(ilwu_nafta, rdd_nafta, wvs_nafta)

par(mfrow=c(1,1))

barplot(naftatable, beside=T, main="Respondents' 'Confidence' in NAFTA",ylim=c(0, 0.7),
        xlab=c(""), xaxt="n", col=c("black", "dark grey", "white"), width=c(0.2, 0.2), cex.main=1, space=c(0,1))
axis(side=1, at=c(0.5, 1.3, 2.1, 2.9), labels=c("A great deal", "Quite a lot", "Not very much", "None at all"),
     tick=FALSE, cex.axis=0.8)
legend("topleft", #pch=15, 
fill=c("black", "dark grey", "white"), border="black", col=c("black", "dark grey", "white"),
       legend=c("ILWU", "RDD", "WVS"), cex=0.8)


##This code produces Table 3
##main NAFTA model
library(Zelig)
library(ZeligChoice)
nafta.fit<-zelig(as.factor(nafta)~ilwu + A + B + white + female +
              as.factor(local) + age + I(age^2) + edu + date + hispanic, data=match.data(m.out1),
              weights = match.data(m.out1)$weights, model="oprobit") 
summary(nafta.fit)

##1=a great deal, 2=quite a lot, 3=not very much, 4=none at all

#Simulations and first difference plots
#This code produces Figure 7 
##RDD vs ILWU 
x.ILWU <- setx(nafta.fit, ilwu=1, A=1, B=0, white=1, female = 0, hispanic=0)  
x.RDD <- setx(nafta.fit, ilwu = 0, A=0, B=0, white=1, female = 0, hispanic=0)
x.cas <- setx(nafta.fit, ilwu = 1, A=0, B=0, white=1, female = 0, hispanic=0)

s.out.1<-sim(nafta.fit, x=x.ILWU, x1=x.RDD) 
s.out.2<-sim(nafta.fit, x=x.ILWU, x1=x.cas) 

plot(density(s.out.1$qi$fd[,4], adjust = 2), xlab="probability", xlim=c(-.1,.45),
     ylim=c(0,15), lwd=1.5,
     main = "Implied effect on NAFTA confidence, ordered probit")
lines(density(s.out.2$qi$fd[,4], adjust = 2), col="blue", lwd=1.5, lty=2)

text(0.21,14.8,"Pr(none | Class-A) - Pr(none | Casual)", cex=.7)
text(.36,11,"Pr(none | Class-A) - Pr(none | RDD)", cex=.7)
abline(v=0, lty=3, col=grey(.3))



######
## matching by age and union tenure; comparing A/B attitutudes toward trade 
## This code produces Figure 3
AB<-subset(comb.data,A==1 | B==1)

m.out2<-matchit(A~age + female + white + hispanic + edu + mode, method="genetic", discard="both", data=AB)
m.out3<-matchit(A~ilwu_yrs + female + white + hispanic + edu + mode, method="genetic", discard="both", data=AB)
summary(m.out2)
summary(m.out3)


#difference by age

#This code produces Figure 5
A2<-subset(match.data(m.out2), A==1)
B2<-subset(match.data(m.out2), B==1)

age_match<-rbind(table(A2$imports_cat)/(nrow(A2)), table(B2$imports_cat)/(nrow(B2)))

par(mfrow=c(1,2))

barplot(age_match, beside=T, main="Support for Import Restrictions \n  Class-A/Class-B Members Matched on Age",
        ylim=c(0, 0.6), xlab=c(""), ylab=c("as % of matched Class-A and Class-B ILWU respondents"), xaxt="n",
        col=c("black", "dark grey"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1)) #, cex.axis=0.6)
#axis(2,cex.axis=0.2)
axis(side=1, at=c(0.35, 1, 1.7), labels=c("Favor", "Oppose", "Don't Know"), tick=FALSE, cex.axis=0.7)
legend("topleft", pch=15, col=c("black", "dark grey"), legend=c("A-member", "B-member"), cex=0.6)

#difference by union tenure
A3<-subset(match.data(m.out3), A==1)
B3<-subset(match.data(m.out3), B==1)

tenure_match<-rbind(table(A3$imports_cat)/(nrow(A3)),table(B3$imports_cat)/(nrow(B3)))

barplot(tenure_match, beside=T, main="Support for Import Restrictions \n Class-A/Class-B Members Matched on Union Tenure",
        ylim=c(0, 0.6), xlab=c(""), xaxt="n", col=c("black", "dark grey"), width=c(0.2, 0.2), cex.main=0.7, space=c(0,1))
axis(side=1, at=c(0.35, 1, 1.7), labels=c("Favor", "Oppose", "Don't Know"), tick=FALSE, cex.axis=0.7)
legend("topleft", pch=15, col=c("black", "dark grey"), legend=c("A-member", "B-member"), cex=0.6)



##Additional specifications for Online Appendix
###including partyID, income, & survey mode as covariates in the structural model
trade.fit.pidinc<-zelig(as.factor(imports_cat)~ilwu + ilwu:A + ilwu:B + white + female + as.factor(local) + age + I(age^2) + edu + date + as.factor(mode) + hispanic + pol_party + income, data=match.data(m.out1), weights = match.data(m.out1)$weights, model="mlogit") 

nafta.fit.pidinc<-zelig(as.factor(nafta)~ilwu + ilwu:A + ilwu:B + white + female + as.factor(local) + age + I(age^2) + edu + date + as.factor(mode) + hispanic + pol_party + income, data=match.data(m.out1), weights = match.data(m.out1)$weights, model="oprobit") 
summary(trade.fit.pidinc)
summary(nafta.fit.pidinc)
  
### Other matching exercises: matching on union family member and using union-only subsample
#####matching on union family member
######Matching
m.out5<-matchit(ilwu~white+as.factor(local)+age+edu+female+ date+hispanic+union_fam, method="genetic", discard="both",data=comb.data, ratio=2)
summary(m.out5)
###This code produces Table A2.1
##imports model
z.out5<-zelig(as.factor(imports_cat)~ilwu+ A + B + white+female+as.factor(local)+age+I(age^2)+edu+date+hispanic+union_fam,
              data=match.data(m.out5), weights = match.data(m.out5)$weights, model="mlogit") 
summary(z.out5)
#with alternate reference category:
alt5<-as.data.frame(match.data(m.out5))
##original coding: 1=favor, 2=oppose, 3=dk (ref cat)
alt5$imports_cat<-recode(alt5$imports_cat, "1=2; 2=3; 3=1") #1=don't know, 2=favor, 3=oppose (ref cat)

z.out6<-zelig(as.factor(imports_cat)~ilwu+ ilwu:A + ilwu:B + white+female+as.factor(local)+age+I(age^2)+edu+date+hispanic+union_fam,
              data=alt5, weights = alt5$weights, model="mlogit") 
summary(z.out6)

#This code produces Table A2.2
#nafta model
z.out7<-zelig(as.factor(nafta)~ilwu + ilwu:A + ilwu:B + white + female + as.factor(local) + age + I(age^2) + edu + date + hispanic + union_fam,
              data=match.data(m.out5), weights = match.data(m.out5)$weights, model="oprobit") 

summary(z.out7)
##1=a great deal, 2=quite a lot, 3=not very much, 4=none at all

#union-only subsample
union.sub<-subset(match.data(m.out1), union==1)

###This code produces Table A2.3
##imports model
z.out8<-zelig(as.factor(imports_cat) ~ ilwu + ilwu:A + ilwu:B + white + female + as.factor(local) + age + I(age^2) + edu + date+ hispanic,
              data=union.sub, weights = union.sub$weights,model="mlogit") 
summary(z.out8)

x1 <- setx(z.out8, ilwu=1, A=1, B=0, white=1, female = 0, hispanic = 0)  
xx <- setx(z.out8, ilwu=0, A=0, B=0, white=1, female = 0, hispanic = 0)

s.out.fd<-sim(z.out8, x=xx, x1=x1)
summary(s.out.fd) #getting risk ratios for paper

#with alternate reference category:
union.sub.alt<-union.sub
##original coding: 1=favor, 2=oppose, 3=dk (ref cat)
union.sub.alt$imports_cat<-recode(union.sub.alt$imports_cat, "1=2; 2=3; 3=1") #1=don't know, 2=favor, 3=oppose (ref cat)

z.out9<-zelig(as.factor(imports_cat) ~ ilwu+ ilwu:A + ilwu:B + white+female+as.factor(local)+ age+I(age^2) + edu + date +hispanic,
              data=union.sub.alt, weights = union.sub.alt$weights, model="mlogit") 
summary(z.out9)

###This code produces Table A2.4
#nafta model
z.out10<-zelig(as.factor(nafta)~ ilwu+ ilwu:A + ilwu:B + white + female + as.factor(local) +
               age + I(age^2)+edu+date+ hispanic,
               data=union.sub.alt, weights = union.sub.alt$weights, model="oprobit") 
summary(z.out10)

##1=a great deal, 2=quite a lot, 3=not very much, 4=none at all
#### END IO-related models

